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We study competition between two biological species advected by a compressible velocity field. 
Individuals are treated as discrete Lagrangian particles that reproduce or die in a density-dependent 
fashion. In the absence of a velocity field and fitness advantage, number fluctuations lead to a 
coarsening dynamics typical of the stochastic Fisher equation. We then study three examples of 
compressible advecting fields: a shell model of turbulence, a sinusoidal velocity field and a linear 
, velocity sink. In all cases, advection leads to a striking drop in the fixation time, as well as a large 

. reduction in the global carrying capacity. Despite localization on convergence zones, one species 

goes extinct much more rapidly than in well-mixed populations. For a weak harmonic potential, 
one finds a bimodal distribution of fixation times. The long-lived states in this case are demixed 
configurations with a single boundary, whose location depends on the fitness advantage. 
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Challenging problems arise when spatial migrations of species are combined with population genetics. The popu- 
. | lation dynamics of a single species expanding into new territory was first studied in the pioneering works of Fisher, 
y ■ Kolmogorov, Petrovsky and Piscounov (FKPP) [IH3]- 

Later, Kimura and Weiss studied individual-based counterparts 
of the FKPP equation |4j, revealing the important role of number fluctuations. In particular, stochasticity is inevitable 
at a frontier, where the population size is small and the discrete nature of the individuals becomes essential. Depend- 



ingon the parameter values, fluctuations can produce radical changes with respect to the deterministic predictions 
dJU- If f(x, t) is the population fraction of, say, a mutant species and 1 — f(x, t) that of the wild type, the stochastic 
FKPP equation reads in one dimension Q: 



d t f(x, t) = Ddlf + sf(l - /) + yjD g f(l - f)£(x, t) (1) 

£^ , where D is the spatial diffusion constant, D g is the genetic diffusion constant (inversely proportional to the local 
population size iVj), s is the genetic advantage of the mutant and £ = £(x, £) is a Gaussian noise, delta-correlated in 
time and space that must be interpreted using Ito calculus d, 0] ■ 

However, many species, from the distant past [§J up to the present, have competed in liquid environments, such 
as lakes, rivers and oceans. Interesting new phenomena arise when population dynamics is coupled to hydrodynamic 
flows [9(. For example, satellite observations of chlorophyll concentrations have identified lon g li ved, segregated 
patches of marine microorganisms off the eastern coast of the southern tip of South America [lfj, where species 
domains are largely determined by the tangential velocity field obtained from satellite altimetry. In cases such as 
these, the dynamics takes place in the presence of advecting flows, some of them at high Reynolds numbers fill ) . 
It is often appropriate to consider a compressible velocity field, both because of inertial effects [l2j associated with 
fairly large microorganisms (diameter 5-500/xm), and because photosynthetic bacteria and plankton often control their 
buoyancy to stay close to the ocean surface [13[. In the latter case, the coarse-grained velocity field advecting the 



microorganisms will contain a compressible component to account for downwcllings [141 ]. Recent works studied the 
dynamics of a single population in the presence of a compressible turbulent velocity field in one and two dimensions 
15, 16]. Here, the interplay between turbulent dynamics and population growth leads to quasilocation on convergence 



zones and a remarkable reduction in the carrying capacity. 

The study of competition in a hydrodynamics context, where both a compressible velocity field and stochasticity 
due to finite population sizes are present, calls for a nontrivial generalization of Eq. (fT]). One complication is that, 
because of compressibility, the sum of the concentrations of the two species is no longer invariant during the dynamics. 
Thus, we must clarify the definition of f{x,t), the fraction of one particular species. A biologically important issue 
arises from overshooting: the density f(x, t) can exceed unity near velocity sinks, resulting in an unphysical imaginary 
noise amplitude in Eq. (fT]). 

In this Letter, we overcome these problems by introducing a new model, designed to explore, for the first time, 
how compressible velocity flows affect the competition between two different species which can die and reproduce via 
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cell division. In absence of an advecting field, the model reproduces the coarsening dynamics of spatial population 
genetics, as shown for one spatial dimension in Fig. la. Of particular importance in population genetics is the 
fixation time Tf, defined as the average time needed for one of the two species to reach extinction. The value of Tf 
is proportional to the total population size in mean field approximation (i.e. assuming uniform spatial mixing) and 
grows to L 2 1 D for a one dimensional system Q. Conversely, Fig. lb shows a typical simulation with our model using 
a high Reynold number synthetic velocity field extracted from a shell model of turbulence [l7j]. In striking contrast 
to the result without advecting flows (Fig. la and Ref. @), one now finds competitions carried out on transient but 
persistent sinks in the velocity field. We find that turbulence not only compresses and reduces the overall population 
density, it also greatly reduces the fixation time Tf. Despite quasi- localization onto velocity sinks, the typical fixation 
time is much shorter than if particles were well mixed. We have also considered a compressible sinusoidal velocity 
field and a confining sink arising from a linear velocity profile. Remarkably, we find that these low Reynolds number 
compressible flows can also dramatically reduce fixation times. In the following, we investigate this phenomenon and 
propose a phenomcnological explanation. 

Eq. ([I} is derived by assuming a fixed total concentration of individuals, so that if f(x,t) is the concentration 
of one species, the other will have a concentration exactly equal to 1 — f(x,t). To describe cases in which the 
total concentration can change due to an advecting flow, we introduce an off-lattice model in which two different 
organisms, A and B, advect and diffuse in space, while undergoing duplication (i.e. cell division) and density- 
dependent annihilation (death). Specifically, we implement the following reactions: each particle of species i = A, B 
duplicates with rate Hi and annihilates with a rate p,ifii, where fij is the number of neighboring particles (of both 
types) in an interaction range 6 oc 1/N, and N is the total number of organisms that can be accomodated in the 
unit interval with total density ca + c B = 1. To reduce the number of parameters, we set pA = Pb = Ms = M; but 
take ha = m(1 + s ) to allow the possibility of a selective advantage (faster reproduction rate) of species A. We will 
start by analyzing in depth the neutral case s = and consider the effect of s > in the end of the Letter. In one 
dimension, our macroscopic coupled equations for the densities CA{x,t) and c B {x,t) of individuals of type A and B 
in an advecting field v(x,t) read 



d t c A =-d x {vcA)+DdlcA+ptc A {l+s-CA-CB)+o-A£, 

d t c B = -d x (vc B )+Dd x c B +HC B (l-c A -c B )+a B ^' (2) 

with a a = \/m c a(1 + s + ca + cbJ/N and a B = \/m c s (1 + c-A + c B )/N. As discussed in [l8j], for an interval of size 
L, let TV -1 — > L/N = S/Ni, where Ni is the local carrying capacity of a region of size S. £(x,t) and £'(x,t) arc 
independent delta-correlated noise sources with an Ito-calculus interpretation as in Eq. ([J) , and we have rescaled the 
particle densities so that the local carrying capacity is 1 when v = s = 0. Although s = may seem nongeneric in the 
context of dynamical systems theory, this is a case of particular interest in population genetics, where it corresponds 
to the neutral theory of Kimura Q . 




FIG. 1: Space-time plot of the off-lattice particle model, (a) no 
velocity field. Simulations are run until fixation (disappearance 
and the much faster fixation time in (b). Parameters values: iV 
in 0. 



advecting velocity field and (b) with a compressible turbulent 
of one of the two species); note the reduced carrying capacity 
= 10 3 , D — 10~ 4 , jU = 1. Parameters of the shell model are as 



181, at this order 



The above equations follow from a microscopic master equation via the Kramers-Moyal method ^ 
equivalent to Van Kampcn inverse system size expansion [19(. Multi-species versions of the FKPP equations have 
been already considered in Rcfs [20j . but not in the presence of an advecting field or number fluctuations. 
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We first simulated the particle model directly for v = s = 0, finding a dynamics similar to that of the stochastic 
FKPP equation, Eq. (JlJ, for s — 0, sec Fig.([TJa) . In this simple limit, our model can be considered as a grandcanonical 
generalization of Eq ([T]), where the total density of individuals ca + cb is now allowed to fluctuate around an average 
value 1. Details of the numerical implementation and the role of density fluctuations are presented in [IsT ]. Hereafter, 
we fix the following parameters: TV = 10 3 , D = 10~ 4 , fi = 1 and L = 1 where L is the one dimensional domain 
endowed with periodic boundary conditions. With these parameters, Tf would be ~ 10 4 for the one dimensional 
FKKP equation, and ~ 10 3 for the well-mixed case. 

Introducing a compressible velocity field v(x,t), as shown in Fig.(^>), leads to radically different dynamics. Indi- 
viduals tend to concentrate at sinks in the velocity field. Further, competition is enhanced and the total number of 
individuals n(t) present at time t is on average significantly smaller than the total carrying capacity TV. We imple- 
mented three different velocity fields: 1) a velocity field v(x, t) generated by a shell model of compressible turbulence 
15 , G3] i reproducing the power spectrum of a high Reynolds number turbulent velocity field of forcing intensity F 
18 1, 2) a static sine wave, v(x) = F sin(27rx), and 3) a linear velocity profile, v(x) = — kx, confining particles near the 
origin. Periodic boundary conditions on the unit interval are used in cases (1) and (2), while in case (3) the system 
size is so large that particles never reach the boundaries. 
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FIG. 2: Average fixation time Tf for neutral competition in compressible turbulence and sine wave advection, as a function of 
(left) the reduced carrying capacity (Z) f and (right) forcing intensity F (smaller (Z) f in the left panel corresponds to larger 
forcing in the right panel). In (a) red circles and blue triangles are particle simulations. Other symbols denote simulations of 
the continuum equations with different resolutions on the unit interval. The black dashed line is the mean field prediction, 
Tf — N(Z)f/2. In (b), only particle simulations are shown and dashed lines are the theoretical prediction Tf — to + c/F based 
on boundary domains, with fitted parameters tq = 9.5, c = 3.5 in the case of the shell model and to = 16, c = 1.4 in the case 
of the sine wave. 



Fig-© shows the average fixation time Tf for s = in the first two cases, at varying the intensity F of advection. 
In the left panel, we plot the fixation times as a function of the time averaged reduced carrying capacity (Z), where 
Z(t) = n(t)/N is the carrying capacity reduction, i.e. the ratio between the actual number of particles and the average 
number of particles TV observed in absence of the velocity field. Plotting vs. (Z)p allows comparisons with the mean 
field prediction, r/ = 2N(Z)p/ fj,, valid for well mixed systems (black dashed line) [f§. For the shell model, we include 
simulations of the macroscopic equations © with different resolutions (256 and 512 lattice sites on the unit interval), 
obtaining always similar results for Tf vs. (Z) p. 

In all cases, the presence of a spatially varying velocity field leads to a dramatic reduction of r/, compared to the 
mean field theory. The fixation time drops abruptly as soon as (Z) < 1, even for very small F. Simulations suggest a 
singular limit of zero intensity (F ~ 0) of the velocity field (and consequently (Z) — > 1 in Fig. 2), as we discuss later. 

These observations can be understood by assuming that global fixation is determined by coalescence of allele 
boundaries like those shown in Fig. (QJ,). Although in the stochastic FKPP with s — the boundaries perform a 
random walk here they are also advected by the velocity field. In particular, the average position of one boundary 
x, neglecting diffusion and number fluctuations, will simply evolve according x = v(x, t). The fixation time is then 
determined by the time needed for the boundaries to reach the center of a sink and annihilate. Because periodic 
boundary conditions are imposed, the number of interfaces is necessarily even, so all of them annihilate pairwise. 
Upon expanding the velocity field around a sink position xq, v(x) = k(xo — x), we find a fixation time ~ k . Hence, 
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we model the average time to fixation r/ vs. the velocity field intensity as Tf=To+ck~ 1 , where To represents the typical 
time to reach fixation at large strain rate k and c is an dimensionlcss constant of order unity. This phcnomenological 
theory describes well the data of Fig. [3J right panel (dot-dashed lines), where in the turbulent case we take the forcing 
F as a proxy for k (k = 2irF for the sine wave and k = J ((Vw) 2 ) cx F for turbulence). 

We now study particles subject to a converging linear velocity, v(x) = —kx. When k is large (and (Z)k is small), 
the fixation time is comparable to the mean field prediction. However, as the forcing decreases, it becomes much 
larger than the mean field prediction (Fig. [3^,, dashed line). In this regime, the fixation time probability distribution 
is bimodal (Fig. [3^, inset): roughly half of the realizations have a short fixation time (faster than mean field), while 
the other half maintain coexistence for an extended period. The space-time evolution of these configurations in Fig. 
[3Jd reveals that the two mutant species are demixed, one on the left and one on the right of the sink. Correspondingly, 
the average total heterozigosity (H(t)), defined as the probability of two random individuals to be of different type 
Q, decays exponentially for very large k (consistent with mean field theory |(|), but tends to a constant, non-zero 
value when smaller values of k are chosen (Fig (3J;) . 
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FIG. 3: (a) Average fixation time t/ in the presence of a linear converging flow v = — kx, as a function of the reduced carrying 
capacity {Z)k- Inset on left shows the distribution of fixation times in the case k — 0.075 and {Z)k ~ 0.38, the right peak 
representing all realizations with fixation times t > 2000. (b) Space-time plot of a realization with a very long fixation time, 
(c) Average heterozigosity as a function of time for different values of k. 

This remarkable behavior arises because of the different boundary conditions for the linear sink, such that the num- 
ber of domain boundaries is no longer always even. Hence, the initial number of boundaries will be odd approximately 
half the time for random initial conditions, leading to the single surviving boundary as shown on the right of Fig. 
This configuration corresponds to a stable stationary solution of the deterministic version of Eqs [5] with s = 0. We 
expect that this demixed solution, with ca{x) and cb{x) nonzero on opposite sides of the sink, should have a lifetime 
(inaccessible in numerical simulations) which grows exponentially with N 0, EH . 




FIG. 4: Coexistence of two species advected by a linear sink k — 0.075, (left) species are neutral as in Fig. 3, (right) Red 
reproduces 30% faster. Notice the shift in the position of the boundary. 
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Finally, what happens for a selective advantage s > 0? Fig. [4] shows a comparison between the neutral case of Fig. 
[3] and a simulation in which red particles reproduce 30% faster (s = 0.3). Even with such a large selective advantage, 
genetic interfaces are still present. In this case, we estimate the shift in the interface as Sx = k~ 1 2\/Ds, by equating 
the outward Fisher genetic wave velocity v g = 2 V Ds [lj with the inward advection velocity at distance Sx from the 
origin. 

To conclude, we have introduced a new model which allows studies of how compressible advection affects Darwinian 
competition. A compressible flow leads to a radically new scenario, with the fixation time being largely determined 
by how species boundaries collapse into the sinks of the velocity field, rather than by diffusive annihilation. Our 
phenomenological theory suggests a significative effect even for very weak compressible fields. Compressible advection 
becomes irrelevant for the fixation time only when the time to drift into the sinks (of order t ~ is much larger 
than the diffusion time (of order L 2 D^, where L is a typical linear size of the population). An obvious application 
of the model is the study of compressible flows in higher space dimensions, including flows relevant to biological 
oceanography of photosynthetic organisms near the water surface. Recent advances in experimental techniques could 
also lead to tests of the one dimensional results. Motility of bacteria has been studied in microstructures of diameter 
comparable or even smaller to that of the organisms j2l| . In this setting, a compressible flow could be created 
by pumping liquid nutrient from the two ends of the tube and extractin g it from the center via a semipermeable 



membrane. Another possibility could be to study floating microorganisms jl3l |22| and make use of the compressible 



effects caused by the vertical component of a converting velocity field [14 1 . 
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SUPPLEMENTARY MATERIALS 



In this note we present the technical derivations of the results in the Letter "Population Genetics in compressible 
flows" ("main text"). First, we describes the derivation of the macroscopic equations for the particle model introduced 
in the main text. Then, we describe details of the continuum simulations and of the shell model of turbulence. Finally, 
we discuss the relation between the model introduced in the main text and the simpler stochastic Fisher-Kolmogorov- 
Pctrovsky-Piskounov (FKPP) equation (Eq. (1) of the main text). 



MACROSCOPIC EQUATIONS FOR THE TWO SPECIES MODEL 

To derive the mean field version of the model discussed in the Letter, we consider first position-independent quantitcs 
n A and n B which represent the number of particles of species A and B respectively. This "zero-dimensional" limit 
corresponds to a single well mixed "deme" or island in a stepping stone model. The birth and death processes are 
characterized by the following transition rates: 



W A (+l,n A ,n B ) = [in A 
W A (-l,n A ,n B ) = fin A (n A + n B ) 
W B (+l,n A ,n B ) = [m B 

W B (-l,n A ,n B ) = fin B (n A + n B ) (3) 

where VF^+L n A, n B ) is the probability per unit time of population A to increase its size by one unit when the total 
populations are n A and n B , and similarly for W A {— 1, n Al n B ) and W B (±l,n A ,n B ). 

There exist several methods to derive macroscopic equations for particle models @, E, HI, most popular being 



the Van Kampen system size expansion and the Kramers-Moyal expansion. While the former is more rigorous 
when studying problems beyond the linear noise approximation, the latter is normally used because it leads to more 
transparent expressions. Moreover, it can be shown (see e.g. [7|, page 251) that at the order of the Fokker Planck 
equation the two methods are completely equivalent. The Kramers-Moyal method works as follows: consider a master 
equation, and call W(An,n) the transition rate from n particles to n + An particles. When n is typically large and 
An relatively small (±1 in the case of birth-death process considered here) a Taylor series expansion of the master 
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equation in An for the probability P(n, t) of having n particles at time t is appropriate: 

d t P(n,t) = [ —^di[ aj (n)P(n,t)} (4) 

with 

aj (n) = J dAn (AnyW(An,n). (5) 

Truncating the Taylor expansion at the second order in j leads to a Fokker-Planck equation. 

It is straightforward to apply the above procedure to two species A and B and obtain the zero-dimensional Fokker- 
Planck equation for the joint probability distribution function P(n A ,n B ,t). Substituting the expression of the tran- 
sition rates of Eq. ([3]) into (@J [S]) leads to: 

d t P{n A ,n Bl t) = -d A [n A ((i - fln A - p,n B )P] - d B [n B {n - P>n A - jln B )P] + 

+ 2^i[ n ^0 u + ^ nA + f ln B)P] + -^d 2 B [n B {^ + fm A + p,n B )P], (6) 
corresponding to the stochastic differential equations: 

dn A {t) 



dt 
dn B (t) 
dt 



n A {n - \m A - \m B ) + \J n A (fj, + p.n A + p.n B )ti{t) 



n B (/j, — fin A — p,n B ) + \Jn B {ii + fln A + fm B )t;'(t) (7) 



where and are delta-correlated Gaussian processes in time with unit amplitude, e.g. (£(^i)£(*2)) = S(ti — t2) 
and are uncorrelated with each other 1 _ (£(ii)£'(i2)) = . For consistency we have to adopt the Ito prescription for 
interpreting the nonlinear noise term [7| . 

From Eq. we see that typical values of n A + n B near the steady state are of order N = fi/jl. Upon defining the 
rescaled particle densities 



n A fi n B fi 

ca = -T7 = -n A , c B = — = -n B , (8) 
N fi N [i 



we obtain 



dc A (t) hic A (l + c A + c B ) 
—jj— = (j,c A (l - c A - c B ) + y — £(t) 



dc B (t) Ulc b (1 + c a +c b ) 

— = fj,c B (l c A c B ) -\- Y ^ C (*)• (9) 

These coupled Langevin equations approximate well the master equation when N = n/p, is large. If N is not large, 
the noise cannot be approximated by a Gaussian noise. 

The model is generalized to include advection and diffusion in space as follows: Particles belonging to the two 
species advect and diffuse in space in a Lagrangian way, i.e. the position a; of a given particle evolves according to: 

d.T , 

v(x) + V2D((t), (10) 



dt 

where f (t) is a unit amplitude Brownian noise source that implements diffusion with diffusion constant D. At each 
time step dt, a given particle duplicates with probability fj.dt or annihilates due to competition with rate p,(h A +n B )dt, 
where h AiB are now the number of particles of type A, B in a region of size S centered in the particle itself. 

The macroscopic equations in the spatial case can be derived using arguments similar to the well mixed case. We 
assume that binary reactions can occur only when both particles are within a length 5, and also take this as the size 
of our discretization volume. Although the total particle number in this region fluctuates, 5 plays a role similar to 
an island spacing in a stepping stone model. We first write the discretized stochastic differential equation for the n A 
and n B , the number of particles in the i th cell of size 6: 



dn A 
~dT 
dn B 
~dT 



n l A (pL — p,n A — fln B ) + diffusion/advection + J n' 4 (u + p,n l A + £in l B )£, A 



n B (n — fin A — fin' B ) + diffusion/advection + W n B (/j, + fin A + (H 
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where the noise terms now obey = S(t — t')Sjj'Sk,k' ■ 

Here, we assume the ration Ni = /i/p is adjusted to match the overall carrying capacity of the i th cell, and that 6 
is big enough to insure the local carrying capacity iVj S> 1. When N is sufficiently large, the diffusion and advection 
terms enter only in the deterministic part of our equations. 

There are two possible choices for passing to the concentration equations: 

1. Change variables to the integrated dimensionless concentration inside each cell: c A B = n A B /Ni. 

2. Change variables to a normalized real concentration (number of particles per unit volume): c l A B = n AB /(N5 d ) 
in d dimensions (we focus here on d = 1). 

In the first case, one is left with a discretized set of well defined Langevin equations, a natural generalization of the 
discrete stepping stone model with dimensionless species concentrations. The first procedure, after setting as before 
Nip, = /i, results in: 



= ^A(l - A - 4) + diffusioa/advection + J *ti±^A 
at V Ni 



^ = ^4(1 -c\- 4) + diffusion/advection + J \f A + &b) & (12) 
at y N[ 

which is a set of well defined stochastic differential equations, valid for a particular lattice spacing S. 
The second procedure, however, leads to a more convenient passage to the continuum limit, 



dc A _ J p.. r.ATVJi J M , ^^U^^ , + 



= c A [fi- pN5(c A - c B )} + diffusion/advection + y — N $ — C 



^ = 4[/z - /W<5(<4 - 4)] + diffusion/advection + ^ ^i±ffl^±M^ (13) 

Starting with Eq. ([13p . a formal way to achieve the continuous limit is to set first pN5 = /x by rescaling the densities. 
This choice corresponds to an average number of particles equal to N in a cell in the absence of fluid flow when the 
total density c A + c B is constant and equal to 1. Then, we take the continuum limit at fixed n, leading to: 

d t c A = ~d x (vc A ) + Dd 2 x c A + /j,c a {1 -c A - cb) + &a£, (14) 
d t c B = -d x (vc B ) + Dd\c B + f.ic B (l -c A - c B ) + <r B £' (15) 

where now < £{x,t)('(x' ,t') >= 0, < £(x,t)£( x' ,t') >= 5( t - t')S(x - x') and < f (x, t)£'(x', t') >= 8(t-t')S(x-x') . 
We also defined a A = ^ fic A (l + c A + cb)/N and a B = \/ fJ,c B (l + c A + c B )/N, We have chosen units of length such 
that L = 1. To return to dimensional units, simply let 1/N — > L/N in a A and a B . To understand the connection 
with stepping stone models, it is helpful to note that L/N = 6/N, so that the amplitude of the noise terms in Eqs. 
(|T4| and (fTB"j) can be written in a similar form as the "generic diffusion constant" discussed in Rcf. Q. Eqs. (fTl|) 
and (|15|) correspond to the macroscopic equations presented in the Letter in absence of a selective advantage, s = 0. 
The derivation of the non-neutral is a straightforward generalization in which the parameter s controls the relative 
difference of the duplication rate of species A with respect to that of species B. Mathematically, this correspondly to 
repeat the expansion but substituting to the first equation in ([3]) the more general expression: 

W A (+l,n A ,n B ) = Wa — > W A (+l,n A ,n B ) = £t(l + s)n A (16) 

leading to the the expression of Eq. (2) in the main text. 

The continuous limit of such reaction-diffusion-advection processes is only a concise summary of the stochastic 
dynamics; underlying this is always the space-discretized version (see Ref. 0, page 314). One reason is that advection 
and diffusion terms can contribute significantly to the variance of the noise when the cells are small (0, page 313); 
these terms can be neglected only if the equations will be discretized on a sufficiently coarse lattice. Moreover, the 
validity of the Kramers-Moyal (or Van Kampen) expansion relies on JVj > 1. Even in absence of a carrying capacity 
reduction due to compressible flows, the number of particles in a sufficiently small cell may not be particularly large. 
The macroscopic continuum equations will then be consistent with results from particle simulations only when there 
is a sufficiently large number of particles when we coarse-grain to scales small compared to the typical scales of the 
density gradients. This is the case for the parameter ranges considered in the main text. 
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CONTINUOUS SIMULATIONS AND SHELL MODELS 



In addition to particle simulations, Eq.s fj L 441 1 5 [) (i.e. Eqs. (2) in the main text) have been simulated numerically by 
using finite difference methods in space and the Euler-Cauchy method for time steps with a variety of velocity fields 
v(x, t). The effect of number fluctuations (i.e. genetic drift) has been taken into account by using the method proposed 
by Doering et. al. [25j : at each time step we compute the position-dependent quantities oa.b = ^ca,b(^ + ca + cb)/Ni. 
If <ja,b > (i.e. if there is a non-zero number of particles of the appropriate type in the region of size <5), then we 
add the noise term; otherwise we set the noise term equal to 0. It has been shown by Doering ct. al. that this 
numerical scheme converges, in a suitable norm, to the continuum limit. In our numerical simulations, we found 
excellent agreement between the results obtained by integrating cq.s (|14M15|) with this technique, compared to those 
obtained by particle methods, as will be discussed in details in a forthcoming paper. 

Next we illustrate how we built the "turbulent" velocity field u(x, t) used in our high Reynolds number simulations. 
Although we focus primarily on one spatial dimension, our analysis of the statistical properties of c(x, t) employs 
velocity fluctuations typical of three dimensional turbulent fluid mechanics. The statistical properties of our advecting 
velocity field u(x,t) will be characterized by intermittency both in space and in time. To obtain a velocity field with 
generic space and time dynamics at high Reynolds number, we build up u{x, t) by appealing to a simplified shell 



model of fluid turbulence 241 ] . Wavenumber space is divided into shells centered around k n = 2™~ 1 fco, n = 1,2, 
For each shell with characteristic wavenumber k n , we describe turbulence by using the complex Fourier- like variable 
u n (t), satisfying the following equations of motion: 

+ Vk*)u n = i[k n+1 U* l+1 U n+2 ~ AfcnWjl-lUn+l 

+ (1 - A)fc„_iU„_iU„_ 2 ] + fn ■ (17) 

The model contains one free parameter, A, and it conserves two quadratic invariants (when the force and the dissipa- 
tion terms are absent) for all values of A. The first is the total energy J2 n \ u n\ 2 and the second is J2 n (~ 1)"^" \ u n\ 2 , 
where a = log 2 (l — A). We fix A = —0.4, as this choice reproduces intermittency features of the real three dimensional 



Navier Stokes equation with surprising good accuracy [241 ] . Using u n , we can build the real one dimensional velocity 
field u(x, t) as 

u(x,t) = Fruity*"* +<(t)e~ tk " x l (18) 

n 

where the constant F controls the strength of velocity fluctuations relative to other parameters in the model. In all 
numerical simulations we use a forcing function /„ = [e(l + i)/u^\8 n ^x, i.e. energy is supplied only to the largest scale 
corresponding to n = 1, where e is the injection rate of turbulent energy. With this choice, the power injected into 
the shell model is simply given by 1/2 J2 n \- U nf n + u nfn\ = e > i- e - ^ i s constant in time. To solve Eqs. (|14j|15j) and 
(|17p we use a finite difference scheme with periodic boundary conditions in space. 



STOCHASTIC FKPP LIMIT 



In this section, we clarify the connection between the advection-free model described in the Letter and in Eqs. 
(|14|15|) and the stochastic FKPP (Fishcr-Kolmogorov-Petrovsky-Piscounov) equation. To handle advecting flows, our 
model gencralizatcs the FKPP equation by relaxing the constraint on the total density, ca(x, t) +CB(x,t) = 1 for all 
values of x and t, see Fig. ([5]). 

To understand this correspondence, it is helpful to change variables from the densities ca, cb to the relative densities 
fA,B(x,t) = CA,B(x,t)/cT{x,t) where cr(x,t) = CA(x,t) +CB(x,t). Eqs. (|14H15j) now become: 

dtfA = -[v-Dd x log(c T )}d x fA + Dd 2 J A + a fA C 

d t fs = -[v - Dd x log(c T )]d x f B + Dd 2 J B + a fB ? (19) 

where now the noise amplitudes are given by a"j A j b = + Ct)/a,b(1 — fA,B)/{Nor) and by definition /^(x,t) + 
f B (x,t) = l. 

Upon substituting ct = 1 and in absence of advection [v = 0), each of the fractions evolves according to the 
stochastic FKPP equations for two neutral alleles. 
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FIG. 5: Reaction plane. The dashed line is the steady state line with selective advantage s = 0, ca + cb = 1, where the model 
is equivalent to the stochastic FKPP equation . The trajectory is a zero-dimensional simulation of Eqs. (|9]) with N = 500 and 
/i = l. 

A further issue is to understand if there is a limit in which the term Dd x log ct can be neglected, so that one 
retrieves Eq. (fT9|) without having to impose any constraint. An estimate follows from studying the fluctuations of 
CT{x,t) = ca{x,€) + Cb(x,€), which obeys a closed equation: 



(20) 



where (£(x, t)^(x't') = 8(x — x')8(t — t'). By defining e(x, t) — cr(x,t) — 1 and assuming e « 1, wc obtain a linear 
stochastic differential equation, 



d t e{x, t) = - M e + Ddle + J ?±£(x, t) 



(21) 



By means of a Fourier transform we can diagonalize the diffusion operator and find that the Fourier modes e(k, i) = 
J dx exp(ikx)e(k, t) are independent Gaussian variables with zero mean and variance <J 2 {k) = 2/j,/[N(fj, + Dk 2 )]. It 
follows that 



(\e(x,t)f 



/' 



2/x f +oc dk 1 1 
TV 7.^ 2?r h + Dk 2 ~ N \ D 



(22) 



The fluctuations in e(x,t) will be negligible when N is large. In addition, it is interesting to study 



(l^logcrl 2 ) ~ <|^e| 2 > = 



+fe — dk 2fiD 2 k 2 
_ k 2ttN (i + Dk 2 



(23) 



where we introduced a large k cutoff (related to the discretization size, |A: maa; | ~ <5 _1 ) in the above integral to avoid 
an ultraviolet divergence. With this premise, the above integral is small if D is small, suggesting that the logarithmic 
terms in Eqs. (|19p can be neglected in this limit. 
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